#include <stdlib.h>
#include <stdio.h>


#define MATACC( m, i, j ) (m[(i)+n*(j)])

void lu_factor( double * m, int n  ) {
	for ( int k = 0; k < n-1; ++k ) {
		// Factor col k
		double pivot = MATACC(m, k, k );
		
		// Substract other lines :
		for ( int i = k+1; i < n; ++ i ) {
			for ( int j = k+1; j < n; ++j ){
				double to_sub = MATACC(m,k,j) * MATACC(m,i,k) / pivot;
				MATACC(m,i,j) -= to_sub;
			}
		}

		// For u vec : divide
		for ( int j = k+1; j < n; ++ j ) {
			MATACC(m,k,j) /= pivot;
		}

		for ( int i = k+1; i < n; ++i ) {
			MATACC(m,i,k) /= pivot;		
		}

		//MATACC(m,k,k) = 1;
	}


}


void fill_matrix( double * m, int n ) {
	// fill with positive rand.	
	for ( int i = 0; i < n; ++ i ) {
		for ( int j = 0; j < n; ++ j ) {
			MATACC(m,i,j) = rand()%4;
			MATACC(m,i,j) = MATACC(m,i,j) > 0 ? 
					MATACC(m,i,j):1;
		}
	}
}


void print_matrix( double *m, int n ) {
	for ( int i = 0; i < n; ++ i ) {
		for ( int j = 0; j < n; ++ j ) {
			printf( "%f ", MATACC( m, i, j ));
		}
		printf( "\n" );
	}

}




int main (int argc, char ** argv) {
	int n = 2;
	if ( argc > 1 )  {
		n = atoi( argv[1]);
	}

	if ( argc > 2 ) {
		srand( atoi( argv[2] ));

	}
	double * m = calloc( n*n, sizeof( double ) );

	fill_matrix( m, n );

	if ( n< 10) 	print_matrix(m, n);
	lu_factor( m, n );
	if ( n<10 ) print_matrix( m, n );


	

}
